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ABSTRACT 

We evaluate the horseshoe drag exerted on a low-mass planet embedded in a 
gaseous disk, assuming the disk's flow in the coorbital region to be adiabatic. We 
restrict this analysis to the case of a planet on a circular orbit, and we assume a 
steady flow in the corotating frame. We also assume that the corotational flow 
upstream of the U-turns is unperturbed, so that we discard saturation effects. 
In addition to the classical expression for the horseshoe drag in barotropic disks, 
which features the vortensity gradient across corotation, we find an additional 
term which scales with the entropy gradient, and whose amplitude depends on 
the perturbed pressure at the stagnation point of the horseshoe separatrices. 
This additional torque is exerted by evanescent waves launched at the horseshoe 
separatrices, as a consequence of an asymmetry of the horseshoe region. It has a 
steep dependence on the potential's softening length, suggesting that the effect 
can be extremely strong in the three dimensional case. We describe the main 
properties of the coorbital region (the production of vortensity during the U- 
turns, the appearance of vorticity sheets at the downstream separatrices, and 
the pressure response), and we give torque expressions suitable to this regime of 
migration. Side results include a weak, negative feed back on migration, due to 
the dependence of the location of the stagnation point on the migration rate, and 
a mild enhancement of the vortensity related torque at large entropy gradient. 

Subject headings: Planetary systems: formation — planetary systems: protoplanetary 
disks — Accretion, accretion disks — Methods: numerical — Hydrodynamics 
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Introduction 



Planetary migration is the process by which forming planets undergo significant 
variations of their semi-major axis, as a result of the tidal interaction with the 
protoplanetary disk. The migration of low mass planets {Mp < 10 — 15 M^), called type I 
migration, has long posed a proble m for scenar i os of planetary formation. Believed to 



be systematically inwards and fast (jWard 



1986 



19971 ). the migration of these planets 



all the way to the central object should be much shorter than the disk's lifetime. This 
raised the issue of why many planetary systems, including our own, harbor planets that 
are not very close to their star, and this triggered a lot of theoretical efforts, during the 
last decade, to exhibit mechanisms that would halt, or even reverse type I migration. 
Long thought to be isothermal (which noticeably simplified the analytical studies as well 
as numerical simulations), the tidal response of the disk w as only recently considered i n 



Morohoshi fc Tanakal fl2003h 



a more realistic manner, by including an energy equation, 
noticed that the disk's response was significantly different in an optically thin disk 
(described in the approximation of the shearing sheet, in which there is no net torque 



between the planet and the disk), while 



Paardekooper fc Mellemal (120061 ) observed in high 



resolution, three dimensional calculations that increasing the opacity of the disk could halt 
or reverse the migration of low-mass planets. This result was subsequently interpreted, 
i n the adiabatic limit, as a result of the advection of entropy within the horseshoe region 



( iBaruteau &: Masset 



2008 



Paardekooper fc Papaloizou 



20081 ). giving rise to a new term of 



the so-called horseshoe drag. This term is of cr ucial importance to t 



migration since it may halt or revert n iigration (IBaruteau &: Masset 



leori e s of planetary 



2008 



2008 



Paardekooper fc Papaloizou 



Kiev fc Crida 



20081 ) . While there is clearly a link between the advection 



of entropy and the existence of an additiona 



scales with the radial gradient of entropy (IBaruteau fc Masset 



term to the horseshoe dr ag, which undoubtedly 



20081 ). there is not yet a 



fully self-consistent treatment of the dynamics of the horseshoe region in an adiabatic disk. 
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Even worse, it is far from clear where the additional contribution originates. So far, it has 
been interpreted as due to overdense and underdense regions which appear at the horseshoe 
U-turns as the result of the advection of entropy in a disk that maintains a pressure 
balance. These regions are strikingly apparent in numerical simulatio ns. They end abrupt ly, 



by a contact discontinuity, at the downstream separatrices (see e.g. 



Baruteau &: Masset 



20081 ). It is easy to realize that the contribution arising from the perturbed density in 



these regions scales with the entropy gradient, and has the correct sign (as inferred from 
numerical simulations). Nevertheless, a simple attempt to perform a full horseshoe drag 
calculation by integrating on the horseshoe streamlines reveals a paradox. The standard 
horseshoe drag calculation consists in integrating over the upstream flow, on each side 
of the planet, the mass flow rate inside of the stream tubes, multiplied by the jump of 
angular mor nenturn expe rienced by the fluid particles that execute a horseshoe U-turn in 



these tubes (IWard 



199ll ). The net result corresponds to the rate of exchange of angular 



momentum between the planet and the gas of the coorbital region. On the upstream 
side of the horseshoe U-turns, the mass flow rate is that of the unperturbed disk, and it 
cannot depend on the equation of state of the gas. Therefore, if the torque due to a given 
stream tube depends on the equation of state, this should be due to the fact that the 
jump of angular momentum during the horseshoe U-turn depends on the equation of state. 
For the sake of definiteness, consider the side of the horseshoe region where we expect 
the material to become underdense, and assume that the disk has a vanishing vortensity 
gradient. In order to decrease the density downstream of the horseshoe U-turns, the flow 
has essentially two possibilities: streamlines can move apart (hence horseshoe U-turns are 
not radially symmetric), or the velocity on the streamlines can differ from the Keplerian 
velocity, and in this case be larger than the Keplerian velocity in the corotating frame, so 
as to lead to an expansion of the material. In other words, a fluid particle originally on 
a streamline at —x = r — Vp [r being its orbital radius and Vp that of the planet) is not 
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mapped to x but to a; + ^, and its velocity in the corotating frame is not V = 2A{x + C,) but 

V = 2A{x + C.) + V (here A = {l/2)rdr^/dr is the first Oort's constant, that quantifies the 
Keplerian shear at the planet's orbit). The perturbed elongation ^ and perturbed velocity 

V adjust so as to create an underdense region downstream of the horseshoe U-turns. A 
simple relationship between and v can be found using the conservation of the Bernoulli 
invariant in the corotating frame. For the sake of simplicity, we assume that here there is 
no pressure gradient, hence the variation of enthalpy during the horseshoe U-turn cancels 
out for an adiabatic flow. The conservation of the Bernoulli constant therefore reduces to 
the conservation of the Jacobi constant J, which reads J = + 2QAx^. To lowest order 
in V and ^, the conservation of the Jacobi constant yields 2BC, + v = (where B = Q + A 
is the second Oort's constant). The jump of angular momentum of fluid elements is, to 
lowest order, Aj = 2rpBx + 2rpB{x + ^) + VpV, hence it is Aj = AvpBx, the very same 
value of the jump that is used in classical, barotropic estimates of the horseshoe drag. 

It seems therefore that there is no room, in a classical horseshoe drag formulation, for a 
dependence on the thermodynamics of the gas, whereas such dependence is evidenced by 
numerical simulationj^. In this work we address this apparent contradiction. We lay out 
our assumptions, and deflne our notation and conventions in section [2l We introduce a 
formalism that leads to a rigorous horseshoe drag expression in section [3l This expression 
features a dependence on the entropy gradient, as expected. We call the corresponding 
term the adiabatic torque excess. We compare this adiabatic torque excess to results of 
numerical simulations in section HI We present the features of the horseshoe region that 
most differ from the barotropic case in section [5l These are essentially the appearance of 
vorticity sheets at the downstream separatrices and a mild production of vortensity all over 

^Although our simple example is limited to a case with no pressure gradient and no 
vortensity gradient, we also expect an adiabatic torque excess in this case. 
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the downstream sides of the horseshoe region. We also work out the pressure and density 
response over the coorbital region, and identify the term that is related to the adiabatic 
torque excess. In section [HI we give a simple interpretation of the origin of the adiabatic 
excess, and we provide suitable torque expressions in section [3 We then discuss additional 
issues in section El and draw our conclusions in section [9l 



2. Prerequisite 
2.1. Notation 

We consider a planet of mass Mp orbiting a star of mass M*, on a fixed circular orbit of 
radius a and of angular frequency Vtp. The planet is immersed in a gaseous protoplanetary 
disk, such that the planetary orbit is coplanar with the disk and prograde. We use S to 
denote the disk's surface density, P to denote the vertically integrated pressure, and T to 
denote the vertically averaged temperature. We assume that the disk's gas follows the ideal 
gas law, which reads: 

p^^. (1) 

/i 

in which TZ is the ideal gas constant and /i is the mean molecular weight. The disk 
isothermal sound speed is therefore: 



cr = w^, (2) 



while adiabatic sound waves propagate at the speed: 



(3) 

where 7 is the usual adiabatic index. We shall consider the gas entropy in two flavors. 
The first of those, which we denote with a lower s, is oftentimes used in the astrophysical 
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literature owing to its simplicity, as its reads: 

Tn ^ ' 

Nevertheless, we shall also need an expression for the entropy that is compliant with the 
original definition of entropy (the entropy exchanged is the heat exchanged divided by the 
temperature). We denote the latter with an upper S: 



Another thermodynamics variable that we shall need is the enthalpy 77, which reads for an 
ideal gas: 



In what follows we shall assume, for the sake of simplicity, that TZ/ ix — 1, which amounts 
to changing the units of temperature, entropy and enthalpy, without loss of generality. 

We note Tq the temperature of the unperturbed disk, which we assume to be uniform. 

We identify a location in the disk by its distance r to the star and its azimuth with 
respect to the planet. The gas has a radial velocity Vr and an azimuthal velocity in 
the frame corotating with the planet (hence its angular frequency in an incrtial frame is 
Vt = v^/r -\- VLp). The radius of corotation Tc is the location in the unperturbed disk where 
the material has same angular frequency as the planet: il(rc) = ilp. We will make use of 
the distance x to corotation: x = r — Tc- We also consider the disk pressure scaleheight 
H = Cg/fl, and the aspect ratio h = H/r. 

The gravitational potential exerted on the disk can be decomposed as the sum of the 
stellar potential = — GM^./r, of the planetary potential $p and of the indirect potential 
The expression of these last two terms are respectively: 





GMp 



(7) 



(H — 2ra cos + + 6^)1/2 ' 
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where e is the softening length of the planetary potential, and 

GM„ 



-r cos = qaVt r cos 



where q = Mp/M^. 

Using the notation of 
across corotation as: 



Baruteau fc Masset 



(l2008l ). we define the gradient of vortensity 



^, _ (ilog(S/cu) _ 3 

dlogr 2 Tc dr 



(9) 



where u = {l/r)dr{r'^Q) is the vertical component of the flow vorticity, and where Sc is the 
unperturbed surface density at the corotation radius. Note that instead of the vortensity 
properly speaking, which is cu/S, we shall oftentimes consider its inverse, S/cu. One reason 
for this choice is that, if we consider that the flow is that of the unperturbed disk, the 
perturbations of the inverse of the vortensity can directly be converted into perturbations 
of density, whose impact on the torque is straightforward. Finally, we define the gradient of 
entropy across corotation as: 



S 



1 d log So 7 — 1 
7 dlogr 7 



rdrSo- 



(10) 



2.2. Basic equations 

The governing equations of the flow are the equation of continuity, the Euler equation 
and the energy equation, together with the closure relationship provided by the equation of 
state. The equation of continuity reads, in the frame corotating with the planet: 

dtJ: + -dri^^rvr) + -d^i^v^) = 0. (11) 

The Euler equations read, respectively in r and 0: 

dfVr + VrdrVr + —O^Vr — T^ll — 2QpV^ = — — 9^$, (12) 
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and 



DtJ 



(13) 



where Dt = dt + Vrdr + and j — r^Jl is the specific angular momentum. The equation 
of state that we adopt is that of ideal gases, which reads P — ET. In that case, the internal 
energy density is e = p/(7 — 1), and the energy equation reads: 



This equation does not include source or sink terms of energy, as we assume the flow to be 
adiabatic. 



Since different authors represent the horsehoe region in different ways, we find it 
useful to refer to properties of this region in a manner independent of the representation. 
Orienting the azimuth rotation- wise, and assuming that, in the corotating frame, the planet 
is located at azimuth = 0, we say that something occurs in front of the planet if it occurs 
at > 0, and behind or on the rear side of the planet if it occurs at < 0. Similarly, we 
refer to the upstream part of the horseshoe streamlines as the set of fluid elements that 
have not performed yet a horseshoe U-turn, whereas the downstream part correspond to 
the set of fluid elements that have already performed their U-turn. 




(14) 



2.3. 



Conventions 



2.4. Assumptions 



Our main assumptions are as follows: 



• We assume that there is only one stagnation point in the vicinity of the planet. 
This is usually not the case in barotropic situations, in which one has generally two 
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X-sta gnation points on the corotation, on each side of the planet (see e.g. iMasset et al. 



2006 



or paper I). In the adiabatic case, nevertheless, we usually have only one 
stagnation point, unless either the entropy gradient is very small, or the potential's 
softening length is very small. We shall therefore assume in what follows that there is 
no ambiguity when we refer to the stagnation point. We will make some additional 
remarks on this assumption in section 18.11 

We assume that the initial temperature field is uniform. This assumption is very 
similar to the assumption of global isothermality that is required in the isothermal 
case to carry out a rigorous horseshoe drag calculation (see paper I). This assumption 
is required in order for the quantity that we shall define in section 12.51 to be conserved 
along a fluid element's path. We will mention in section how our results can be 
generalized to the case of an arbitrary temperature profile. 

Finally, we assume that the flow is in steady state in the frame corotating with the 
planet, and that the flow upstream of the horseshoe U-turns is unperturbed. We note 
that these two assumptions are contradictory: reaching a steady state implies that 
the flow has executed many horseshoe librations, hence the upstream flow cannot be 
that of the unperturbed disk. Nevertheless, we restrict ourselves, in this analysis, to 
a study of the unsaturated torque value. This implies that we consider the flow on a 
time scale (i) longer than what it takes to execute a horseshoe U-turn, so that the 
flow has reached a steady state in some region of interest enclosing the planet, (ii) 
shorter than the horseshoe libration time, so that in the region of interest any fluid 
element has executed at most one horsehoe U-turn. 
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2.5. A useful invariant 

We consider the Jacobi-like Bernoulli constant whose expression is, in steady state: 

Bj = ''"^^-^f+'^r + + + _ 1,2^2 ^ ^_ (15) 

Since Bj is conserved along the path of a fluid element, the quantity defined as: 

G = B.j- ToS (16) 

is also conserved in an isentropic flow. Far from the planet, assuming a purely azimuthal 
motion, we can write, from Eqs. (fTSl) and f|T6l) : 

drG ={n- np)ruj + (T - To)drS, (17) 

where we have used the rotational equilibrium, which can be deduced from Eq. ( fT2l) by 
letting Vr = 0, and which reads: 

dr(j), + drV = m'^ + TdrS. (18) 

As a consequence, in the unperturbed flow (T = Tq), G is maximal at corotation, as is Bj 
in an isothermal disk, and it admits the following expansion to second order in x: 

G = Gc + 2ABx^, (19) 

where Gc is the value of G at corotation, A = {l/2)rdQ/dr is the flrst Oort's constant 
and B = {l/2r)d{r^Q)/dr is the second Oort's constants, which are to be evaluated at 
corotation. 

3. A horseshoe drag expression 

3.1. Initial formulation 

We consider a domain D of the disk consisting of an angular sector centered on the star 
that encloses the planet, as depicted in Fig. [H The azimuth of the boundaries (1) and (III) 
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is chosen sufficiently large so that the material that crosses them can be considered to 
have a purely circular motion, while the boundaries (II) and (IV) enclose the horseshoe 
region, and are chosen much further away from corotation than the horseshoe separatrices, 
so that they can be considered as circular. We note that in a real situation, this last 
requirement is hindered by the wake, which corresponds to epicyclic motion excited at 
Lindblad resonances. Nevertheless we discard this behavior in the present work, which 
amounts to considering an isolated corotation region. We want to avoid flux of material 
into the domain D through the boundaries (II) and (IV), so we adopt streamlines to define 
them. These two boundaries are therefore associated to some value of the G invariant. For 
convenience, we use the same value for G on both sides, and we note it GocC. 

The torque exerted on the planet by the material included in the domain D reads: 

^ = j j ^id^^)d(j)rdr, (20) 

We can transform the integrand of Eq. fl20l) as follows: 

E(9^<l>)r = -rd^P - d^iv^T.]) - dr{rT.Vr3), (21) 

where we have made use of the assumption of a steady state in the corotating frame, and 
where we have used Eqs. (ITT]) and (IT^ . 

Since the boundaries of the domain D are sufficiently far from the planet that we can 

^ Since G = 2ABx^ to lowest order in x, it is certainly possible to adopt the same value 
of the G invariant for the outer and inner boundaries of the integration domain, if these are 
not too far from corotation. If they are located at a sizable fraction of the semi-major axis 
from corotation, however, this may not longer be possible, and one may have different values 
of the G-invariant for the inner and outer boundaries. This does not change the generality 
of the demonstration presented here. 



Fig. 1. — Sketch of the domain of integration (hght grey) enclosing the horseshoe region 
(dark grey). The dashed hne shows the corotation. As can be seen here, there is usually 
an offset between the planet, and the stagnation point of the horseshoe region (the point of 
intersection of the boundaries of the horseshoe region, i.e. the separatrices) . There is also 
generally an offset between the corotation radius and the orbital radius of the planet, but it 
is much smaller. 
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neglect radial motions, Eq. fl20l) can be recast as: 



[rP + r{n-np)Ej]dr 



R 



(22) 



where (r^) is the radius of the outer (inner) streamhne that has G = Goo- The R and 
F notation at the right bracket of Eq. fl22l) means that the integration over r has to be 
performed respectively at the rear or front side of the domain D (boundaries III and I in 
Fig. [T]). In Eq. ( l22l) . the first term of the integrand represents the pressure torque exerted 
on the material enclosed within the domain D, while the second term represents the budget 
of angular momentum brought to this region by advection. Since the flow is steady in 
the corotating frame, the angular momentum of this domain is constant in time and the 
torque is therefore integrally transmitted to the planet. We recogniz e in the second ter m 



of the integrand of Eq. ((221) the classica l horseshoe drag expression ( Ward 



Masset 



2001 



Masset fc Papaloizou 



2003 



Masset et al. 



2006 



1991 



1992 



Paardekooper fc Papaloizou 



2009al ). but this equation also shows the pressure contribution, which has been overlooked 



in previous analysis. 



3.2. A simplifying assumption 

We assume that in the unperturbed flow, far from the planet, lines of constant entropy 
and of constant G coincide. This is not exactly true: in the unperturbed disk the lines 
of constant entropy are circles centered on the primary, while, owing to the presence of 
the potential's indirect term given by Eq. ([8]), the iso-G lines have a variable distance to 
corotation in the horseshoe region (an extreme example being provided by the tadpole 
separatrix). Our assumption is therefore more appropriate of a situation in which we 
discard the indirect term of the potential, or of a shearing sheet situation. Nevertheless, 
we shall make hereafter this assumption, anticipating that our torque expression does not 
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depend sensitively on it. 

This assumption allows to define unambiguously the derivative dcS, which can be 
obtained from the measure of S and G on two neighboring streamlines, since these variables 
are intrinsic to the line. The quantity dcS is therefore also intrinsic to the streamline, and 
hence conserved along it. 

In the unperturbed disk, this derivative has the following expression: 

It therefore diverges at corotation. 

Mass conservation yields a relationship which will be useful to evaluate the production 
of vortensity during horseshoe U-turns. We consider a horseshoe stream tube, of width 6G. 
The mass flux across the tube is: 

SG 

M = j^r{n - (24) 
We can transform Eq. (JT71) to write: 

{1-6T- dGS)drG = {yi- Vtp)ru, (25) 
where 6T = T — Tq. This eventually gives the mass flux: 

M = 5G-{l-5TdGS). (26) 

OJ 

We can therefore write the following relationship: 



UJ 



;i - 5TdGS) 



d ^ 



(27) 



where the subscript d {u) respectively denote the quantities downstream (upstream) of a 
horseshoe U-turn. 
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3.3. An alternate horseshoe expression 



We now perform a change of variable in the integrals of Eq. fl22l) . We choose as the 
new variable of integration the value of the G invariant. We have therefore to split the 
intervals of integration into intervals over which the dependency of G on r is continuous 
and monotonic. Anticipating on the results exposed in the following sections, we note that 
the value of the G invariant is not necessarily the same for the front and rear separatrices, 
as shown in Fig. [2J We denote their respective values with G^ and G_. 

Using Eq. (I25l) . Eq. fl22|) is therefore recast as: 



We split the integration on the downstream side of the horseshoe U-turns in order to 
account for the possible discontinuity of G at the separatrices. More precisely, the first 
integral over G in Eq. (128!) corresponds to an integral over intervals (1) and (2) (see Fig. [2]), 
the two following integrals (from Gc to G_ and G^ to G^o) correspond respectively to 
intervals (3) and (4). In a similar manner, on the front side, the intervals of integration 
over G are, in order of appearance, (5), (6) and eventually (7) and (8) together. We also 
note that in Eq. fl28l) we have substituted the total pressure P by the perturbed pressure p, 
since the unperturbed pressure field is axisymmetric and does not contribute to the pressure 
torque. 



r 




F 



(28) 
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> 

Azimuth 



Fig. 2. — Sketch of the flow topology that shows the value of the G invariant along the 
separatrices. The value of the G constant is continuous across the upstream separatrices, 
while it displays a jump across the downstream separatrices. The different intervals of 
integration mentioned in the text are shown along the y-axis. Intervals (1) and (5) are 
meant to extend up to r^, and intervals (4) and (8) are meant to extend up to r^. 
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Other comments are in order regarding the transformation of Eq. fl22|) into Eq. fl28l) . 
The integrand in Eq. fl22]) is regular, and ahhough it may display discontinuities (we shall 
see that Q is discontinuous at the downstream separatrices, as well as S) it does not 
contain any singularity. In contrast, the integrands of Eq. (128|) feature the vorticity, which 
is singular wherever Q is discontinuous. Since the dependency of G upon r on the different 
intervals of integration is regular, the integrals on G must exclude any singularity (simply 
because the contribution of the integrand of Eq. fl22l) on any vanishingly small interval in 
r is vanishingly small, and the same must be true of the equivalent integral in G). As a 
consequence, any integral of Eq. ( I28i) which has a boundary corresponding to a downstream 
separatrix must be understood as excluding the possible singularities at this boundary. For 
instance, the second integral over G in Eq. (l28l) must be understood as: 

/ (l~STdGS)-jdG= lim / (l-STdGS)-jdG. (29) 
Jgc ^ SG^o+Ja^ uj 

We also note that including such singularities would be mathematically ill-defined anyway, 
because they lie on the edge of the integration domain. 

Eq. (1281) can be further simplified as explained below. The specific angular momentum 
j of a fluid element that has some value of the G invariant can be written as: 

J = Jo + 5 J, (30) 

where Jq is the specific angular momentum of fluid elements in the axisymmetric, 
unperturbed flow which have same value of G, and 6j is the difference, that we wish to 
evaluate. Azimuthally far from the planet, Eq. f|T6|) can be recast as: 

G = + $^ _ + ^ _ ToS, (31) 

where we have neglected the indirect term (as discussed in section 13121) . Differentiating 
Eq. fl3T|) . we can write, to first order in the perturbation: 

= r6r{n^ - 2mp) + r^5n{n - fip) + dr^Jr + (32) 
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The left hand side of Eq. fl32|) cancels out because, as stated above, we consider a fluid 
element that has same value of G in the unperturbed flow and in the perturbed one. Using 
the rotational equilibrium of the unperturbed flow, and writing 6P = drPoSr + p, we can 
recast Eq. (132|) as: 

= (fi - np)6j + ^. (33) 

We therefore have an expression for 6j that involves the perturbation of pressure p. 
Since we consider the variation of angular momentum of a fluid element for a constant value 
of G, it is natural to change to G the variable of the integrals of the pressure torque. The 
generic integrand of Eq. fl28p therefore becomes: 

^{l-5TdGS)j + ^ = ^[j:r{n-n,)UG) + J:r{n-%)5j+rp] (34) 

= -a - STOgSMG) + j^^^rp, 

where we have used Eqs. fl25|) and fl33l) . 

The last term of the right hand side of Eq. flM|) cancels out to first order in the 
perturbation. Furthermore, the parity in (j) of the perturbed density and of the perturbed 
pressure (which will be studied in section 15.31) is such that the rear and front contributions 
of this term to Eq. (128!) cancel out. Therefore we can discard the pressure torque in Eq. (128|) . 
provided we use the unperturbed value of the specific angular momentum jo- Since G has a 
maximum at corotation, the notation jo{G) contains an ambiguity, as one has to specify if 
jo has to be evaluated inside or outside of corotation. In what follows we use the notation 
Jq to remove this ambiguity, and we drop the G dependency for the sake of brevity. 

The fact that one can neglect the pressure torque and use the value of the unperturbed 
angular momentum has been demonstrated in paper I for isothermal disks. It is linked 
(i) to the fact that the evanescent pressure waves excited by the perturbation of the fiow on 
the downstream side of the horseshoe U-turns do not alter the total horseshoe drag, even 
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though they redistribute the perturbation of density and therefore change the local torque 
density, and (ii) to the fact that the impact on the torque of the rear-front asymmetry due 
to the feed back of the evanescent waves cancels out. 

Using the simplification described above, we transform Eq. fl28l) as follows: 



G_ 



Jo dG + 
hdG 



G 
Gc 



- s 






0- 


pGc 


S 


G+ 


UJ 



J^dG + 
hdG + 



pGoo 


s 








jUg 


Ig+ 




0+ 



(35) 



R 



0+ 



/•Goc 






j+dG 


Jg, ^ 


0+ 



Apart from discarding the pressure contribution and changing j into jo in Eq. fl35l) . we 
have introduced an additional simplification by making use of Eq. f l27|) . and by referring 
to the unperturbed, upstream value of T,/uj{l — STdcS). We explicitly specify by the use 
of the index 0^ whether the unperturbed vortensity must be considered outside or inside 
corotation. 

We note (see Fig. [2]) that the integrals over the intervals (1) and (5), on the one hand, 
and over the intervals (4) and (8) on the other hand, exactly cancel each other, hence the 
torque expression can further be reduced to: 
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(36) 



where Ajo = Jo^ — Jo • "^^^ corotation torque has therefore a very sim ilar expression to the 



usua 



(barotropic) case, when integrated over a Bernoulli-like variable (IMasset fc Papaloizou 



20031 ). except that the domain of integration does not have the same limits on the front 



side and on the rear side. Denoting by Gg the arithmetic mean of the neighboring values 
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G- and 6*+, we can transform Eq. fl36|) into: 

r 















0+ 





(37) 



where Ajq is the value of Ajq for G = G^, and where AG = G_|_ — G_. The left term of the 
R.H.S. of Eq. fl?r|) is the standard horseshoe drag, which scales with the vortensity gradient, 
and the right term, which we hereafter denote with Fi, is an additional contribution which 
arises from relaxing the barotropic hypothesis. Evaluating this contribution therefore 
amounts to evaluating AG. This is the purpose of the following section. 



3.4. Discontinuities at the stagnation point 

The discontinuity in G that we wish to evaluate arises from the presence, at the 
downstream separatrices, of a contact discontinuity (see Fig. [2]). Formally, an inspection 
of Eq. f|T5|) and f|T6l) reveals that this discontinuity can come, far from the planet, from a 
discontinuity in the enthalpy, in the entropy, and in the azimuthal velocity (the other terms 
being continuous across the separatrices). The situation becomes simpler at the stagnation 
point, where the velocity vanishes and where the discontinuity in G is entirely accounted 
for by a discontinuity of entropy and enthalpy. We have therefore: 

AG = r/+-r/_-To(5+-5_), (38) 

where the + (— ) sign refers to quantities outside, or in front of (inside, or behind) the 
stagnation point. Eq. (135]) can be transformed into: 

AG = (4^^ - s^^)^— Pi^-i)/^ - ^°-(logs+ - logs.), (39) 
7 — 1 7 — 1 

where Pg is the value of the pressure at the stagnation point (this value is well-defined 
at this point, on the contrary to the other thermodynamics variables, since the pressure 
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field is continuous). Using the first order expansions: s\^'^ — s^I'^ (s+ — sS)s\l'^ ^ and 
logs+ — logs_ ~ (s+ — s_)/so, we are led to: 



AG'= - 

7 



l_3l/7-l(p(7-l)/7 _ p(^-^)/^)(,^ - ,_), (40) 



where Pq and sq are the pressure and the entropy at the location of the stagnation point, in 
the unperturbed disk. 

In the particular case of a small perturbation: \Ps — Pq\ <^ Pq^ the discontinuity in G 
takes the form: 

= (41) 

7^0 So 

where ps = Ps — Po is the pressure perturbation at the stagnation point. 

In order to proceed in the evaluation of AG, we use the conservation of the G invariant 
along the upstream separatrices. At this stage of the derivation, we note that numerical 
simulations reveal that the horseshoe region in the adiabatic regime is rather asymmetric, 
so that the distance of the separatrices depend on the quadrant under consideration (rear 
vs. front, and upstream vs. downstream). Denoting with the distance of the upstream 
separatrix from corotation on the rear side, we have, by evaluating G respectively on that 
separatrix far from the planet, by means of Eq. (|T9l) . and at the stagnation point: 

G_ = a + 2AB{x^f, (42) 

and 

G_ = G,-r]c + V- + % + To{S, - S.), (43) 

where $p is the planetary potential at the stagnation point, rjc and Sc are respectively 
the enthalpy and entropy at corotation in the unperturbed disk, and where we have used 
the relationship G^ = $*(rc) - (l/2)r2fi2 + r/^ - ToS^ to derive Eq. From Eqs. 

and fH3|) we deduce: 

2AB{x^y = V--Vc + n + To{S, - (44) 



- 23 - 



In a similar manner, denoting with the distance of the upstream separatrix from 
corotation on the front side, we have: 



2AB{xi )^ = - r/, + + To(5, - S, 



(45) 



From Eq. (jlPj) . we have: 



AG 



7 



So 



1/7 



v± - Vc 



So 



hence we have: 



2AG = -iv^ + r].-2r],y^^^., 

7 So 



(46) 



(47) 



where we have discarded the additional term that scales with S^. Similarly, adding Eqs. 
and (H5|) yields: 

2AB[{x^ f + [x^f] = + r/_ - 2r], + 2%, (48) 

where we have omitted the term Tq{2Sc — S- — 5*+), which scales as K^Vt^x"^, and is hence 
negligible compared to the left hand side. Using Eqs. P7|) and (HHj) . we obtain: 



s+ — s_ 



AG = -{Ai?[(xf)2 + (xr)V<^';} 



(49) 



Since the entropy is conserved during the advection along the upstream separatrices 
we have: 



■50 



F I _R 



The discontinuity of G reads finally: 

AG={AB[{x^f + {x:f]-%}S 
and the associated torque excess is, using Eq. (1371) : 



F I J? 



(50) 



(51) 



Ti = -{xl + x«)^So5 {AS[(x«)^ + [x^f] - %] , 



(52) 
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where we have used the equahty Ajq = 2Ba{x^ + x^). Since Ps > Pq, AG has same sign 



as S, and the torque excess scales with the negative o f the entropy gradient, as a 



observed in previous works (iBaruteau fc Masset 



2008 



Paardekooper &: Papaloizou 



ready 



20081 ) 



We can render Eq. (!52|) a bit more compact by defining the averaged upstream half-width: 

(53) 



Xq 



.X. 



R\2 



[X. 



F\2 



Considering that the relative asymmetry between the rear and front part of the horseshoe, 
albeit larger than typically observed in barotropic situations (see paper I), remains 
moderate, we have: 



Xs-^ix^ + x^), (54) 

and therefore: 

Ti ^ -4x2So5 {2ABxl - , (55) 

This expression displays explicitly the dependency of the torque excess on the entropy 
gradient. As we shall see hereafter, it also brings in a dependency on the potential's 
softening length, as both Xg and $p depend on this quantity. 



4. Comparison to numerical simulations 

We have undertaken a number of numerical simulations in order to test the validity 
of Eq. (1551) . which quantifies the adiabatic torque excess. We briefiy describe the code 
used and the parameters adopted, then we assess the correctness of Eq. (1551) from series of 
calculations in which we vary the entropy gradient or the softening length. 
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4.1. Numerical set up 



We used the code FARGO^I flMassetll2000al ]bl: iBaruteau fc Massetl 120081 ) in its adiabatic 
and isothermal versions. Unless otherwise stated, the polar mesh extends from r^[^ = 0.4 to 
i^max = 1-8, and has a resolution of 1200 x 1200, the zone boundaries being evenly spaced. 
The planet is located at r = 1, and it is held on a fixed circular orbit. Its mass ratio to the 
primary is g = 3 ■ 10^^ {i.e. the planet has one Earth-mass, if the central object has a solar 



mass). The disk aspect ratio at r = 1 is /i = 0.0 5. The adiabatic index o: 



Damping boundary conditions (as described in 



de Val-Borro et al 



the gas is 7 = 1.4. 



20061 ) are used at the 



he frame corotates with the planet, and the Coriolis 



19981 ) . The surface density is a power law of the 



inner and outer edges of the mesh 
force is conservatively implemented (iKley 
radius: S(r) oc r~'^. 

In our first set of calculations, the slope of the surface density is set to cr = +3/2 
(hence the vortensity related part of the corotation torque cancels out), and the potential 
softening length ranges from e = 0.1 to 10, where e = e/H is a dimensionless value of the 
softening length. The logarithm of e is evenly spaced, and 21 calculations are performed, so 
that e = 10*/^*^^^, with < z < 20. These calculations were performed both for an adiabatic 
flow (for brevity we denote them with EAj o<j<2o) and for an isothermal flow (we denote 
them with Elj o<i<2o)- Since the initial temperature field is fiat, the index of the entropy 
gradient is S = cr(7 — l)/7 = +3/7. We therefore expect a negative torque excess in the 
adiabatic case. 

In our second set of calculations, the slope of the surface density ranges from a = —4 
to cr = 4, the softening length being kept fixed at e = 0.3. We have divided the slope 
interval into 41 calculations such that cxj = —4 + 0.2 x i [0 < i < 40). These calculations 



^See: |http : //f argo . in2p3 . f r 
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were performed both for an adiabatic flow (we denote them with SAj^o<j<4o) and for an 
isothermal flow (we denote them with SIj^o<j<4o)- Since the initial temperature field is flat, 
the index of the entropy gradient is 

S = ^^a. (56) 
7 

It therefore ranges from —8/7 to +8/7. 



4.2. A note on the half width of the horseshoe region 

As underlined in section I3.4[ the horseshoe region in adiabatic calculations is more 
asymmetric than in the isothermal case. It also has a different dependence, more complex, 
on the disk parameters. This can be seen in Fig. [3], in which we plot the different values of 
Xs inferred from a streamline analysis on the runs of the series SAj. In particular, the naive 
expectation that 

xf = xT/j'^' (57) 

is true only in the barotropic case (5 = 0). For non- vanishing entropy gradients, the 
averaged upstream half width exhibits a V-shaped dependence on S, whereas the isothermal 
case is insensitive to this parameter. This more complex dependence of the width of the 
horseshoe region has a consequence on the evaluation of the torque excess. A simple 
estimate of the latter is: 

r: = r^ , (58) 

7 

but this estimate assumes Xs^adi = Xg^iso/l^^^, so that one gets rid of the vortensity related 
part of the horseshoe drag, retaining only the excess. This is not exactly true, since the 
adiabatic case has a horseshoe region in general larger than expected. In order to fully 
eliminate the vortensity part of the estimate, one should therefore use the expression: 

Tl = T^-^- ^Vl^JSo [xU - (59) 
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Fig. 3. — Half width Xs of the horseshoe region in the adiabatic case (diamonds), as a 
function of the entropy gradient S. We use the arithmetic mean of and x^ as an estimate 
of Xg. The dotted hne shows x^, and the dashed hne shows x^. The triangles show the half 
width of the horseshoe region (divided by 7^''^) obtained from runs with same parameters, 
but with an isothermal equation of state. The rear/front averaged upstream half width has 
a characteristic V-shape. The data presented here are obtained from the series of runs SA 
and SI. 
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We shall hereafter refer to the estimate given by Eq. fl58l) as the rough or global estimate 
(hence the r superscript), and to the estimate of Eq. fl5^ as the corrected estimate (hence 
the c superscript). 



4.3. Dependence on the softening length 

We compare the torque excess, when e varies, as provided by Eq. ( !55l) and as provided 
by a direct comparison of the torque value in adiabatic and isothermal calculations. The 
results are presented in Fig. HI We see that there is a good overall agreement between the 
direct measure of and the result of Eq. (!55|) . as the relative difference is at most ~ 10 % 
over the range 0.1 < e < 1. The value of Xs, as measured on the output at t = 40 orbits and 
= ±1 rad, is used in the evaluation of Eq. fl55l) . Similarly, the location of the stagnation 
point is determined, and used to evaluate the planetary potential that features in the 
bracket of Eq. fISS]) . We note that in all the runs of this series, there is a unique stagnation 
point in the planet's vicinity, hence there is no ambiguity in its determination. 

There is an approximate scaling of the torque excess in the range 0.1 < e < 1 



(jBaruteau 



20081 ) . Fig. m shows that this trend results partly from the increase of Xs as e 
decreases (dotted curve), and partly from the increase (in absolute value) of the planetary 
potential as the softening length decreases (dot-dashed line). 



4.4. Dependence on the entropy gradient 

We now evaluate how well Eq. fl55|) accounts for the torque excess in the second series 
of runs (SA and SI, see section HTTj) . in which we vary the entropy gradient. 

The results are presented in Fig. O They show a satisfactory agreement between the 
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0.1 1.0 

e 

Fig. 4. — Torque excess as a function of the softening length. The thick dashed hne shows 
the excess directly measured by comparing the runs EA, and Elj, by means of the Eq. (|59l) 
- or fl58l) . since the vortensity gradient is zero in this series of runs while the thick curve 
with diamonds shows the prediction of Eq. (!55l) . The dash-dotted and dotted curves show the 
result of Eq. ( l55l) in which one respectively adopts a constant value for Xg or for $p (namely 
the value measured for e = 1 ) . The grey dashed line shows the prediction of Eq. (1961) , which 
will be derived in section [71 
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Fig. 5. — Torque excess as a function of the entropy gradient. The thick dashed hne shows 
the corrected direct estimate T\ given by Eq. ( l59i) . while the thin dotted hne shows the 
rough estimate T\ of Eq. (!58l) . The thick sohd hne shows the result of Eq. (!55|) . in which 
we use the planetary potential at the stagnation point and the horseshoe width, given by 
a prior streamline analysis performed at t = 40 orbits for each calculation (and corrected 
as explained in the text). The inset plot shows the azimuth of the stagnation point, as 
a function of the entropy gradient (lower axis) and as a function of the corotation radius 
(upper axis). The vertical dotted lines show the zone boundaries. The thick solid lines of 
the inset plot shows linear regression fits of the stagnation point location, aimed at getting 
rid of the spurious oscillations. These fits were used in Eq. (l55l) and gave the thick solid 
curve of the main plot. The grey dashed line shows the prediction of Eq. (|96i) . which will be 
derived in section [71 
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corrected estimate of the torque excess (given by Eq. [59], and represented by the thick 
dashed hne), and the theoretical prediction of Eq. (jSSD, represented by the thick sohd hne. 
This figure also shows that the azimuth of the stagnation point oscillates as the entropy 
gradient varies. The oscillations are significantly large, and can be unambiguously identified 
as a finite resolution issue, for the following reasons: 

• when the entropy gradient varies, the corotation radius is slightly shifted. If one 
displays the azimuth of the stagnation point as a function of the corotation radius, 
the period of the oscillations is equal to the radial resolution. 

• We have carried out higher resolution calculations (2000 x 2000, not shown here, as 
they were performed over a small number of orbits and a steady state was not reached 
in the vicinity of the planet), which were ran for 10 orbits only (the location of the 
stagnation point hardly varies afterwards). They also yield oscillations of the azimuth 
of the stagnation point with a shorter period and smaller amplitude. 

Given the large amplitude of the oscillations of the stagnation point for the resolution that 
we were able to use for this series of runs, we need to smooth the (psiS) function. This 
is done simply by performing a linear regression fit, independently for the rear and front 
stagnation points, as shown in the inset plot of Fig. [51 The fitted value are then used in 
Eq. fl55|) . which yields the solid thick curve. 

The value of 5 = +3/7 used in the series of runs EA and EI corresponds to a location 
where the error due to the oscillations of the azimuth of the stagnation point is small. We 
therefore understand from this plot why there is a good agreement between the measured 
and predicted adiabatic excess in Fig. [H The fact that there is a good agreement for all 
values of e in Fig. [H suggests that the error on the location of the stagnation point is 
essentially due to the location of the corotation with respect to the mesh (which does not 
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depend on the softening length). We recover the fact that, for e = 0.3, |ri|~1.5-10 

It is noteworthy that, on the left part of the plot of Fig. \5[ there is a sizable difference 
between the rough and corrected estimates of the excess (thick dashed curve and dotted 
curve). The reason for this is straightforward, since for these values of the entropy gradient, 
the surface density increases steeply outwards, hence the vortensity related corotation 
torque is large and positive. For large values of \S\, the adiabatic horseshoe half width is 
significantly larger than xf"/'y^^^ (see section l¥!2l) . therefore a fair fraction of the adiabatic 
excess is actually due to the boost of the vortensity related torque. 



We describe in this section the properties of the coorbital region that differ from the 
barotropic case. The most salient feature is the appearance of vorticity sheets along the 
downstream separatrices. We also derive the production of vortensity within the horseshoe 
region. The latter is weak, however, and does not have any impact on the horseshoe drag. 
We finally discuss the properties of the pressure disturbances arising from the horseshoe 
dynamics, and identify which component of these is responsible for the adiabatic torque 
excess. 



The discontinuity of the G-invariant that we evaluated in section 13.41 at the stagnation 
point also exists along the downstream separatrices, as can be seen in Fig. |2j Far from the 
planet, one can write the expression of the G discontinuity across a downstream separatrix 
as: 



5. Properties of the coorbital region 



5.1. Vorticity sheets 




(60) 
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where we have performed the same series of transformations as in Eqs. fl38|) to fHTl) . where 
p'g is the perturbed pressure at the separatrix, and where ^2+ (^^-) is the angular velocity 
on the outside (inside) of the separatrix. If we assume in a first approximation that the 
perturbation of pressure at the separatrix, far from the planet, is much smaller than at the 
stagnation point, most of the G- discontinuity is ensured by the jump of azimuthal velocity, 
i.e. by a vorticity sheet. Denoting with Af = r{VLj^ — Vt^) the magnitude of the vorticity 
sheet (i.e. the singular part of the flow's vorticity reads: uj = Av6{x — Xg), where 6 is Dirac's 
delta function), and assuming that the jump is sufficiently small so that we can write: 

AAx 

+ 2Vip = — , (61) 

Ct 

we are led to: 

A ^ f — , , / „^ at the rear separatrix 
AG I 4„|^|^E. 

; — — - — at the front separatrix 

I 4a|A|a;2So 

Therefore, if there is a positive torque excess (hence S < 0), then there is a negative 
(positive) vorticity sheet in front of the planet (at the rear of the planet), and vice- versa. 
In order to illustrate this, we show in Fig. [6] the vorticity field in a run with 5 = 1/3 
and similar characteristics as those of section 14.11 except that the resolution was set to a 
2000 zones in azimuth and 8000 zones in radius. We measure an adiabatic torque excess 
Fi = —1.14 ■ 10^^^ M^QpTp, by performing a run with similar characteristics and an 
isothermal equation of state. This allows us to predict the jump of angular velocity at the 
downstream separatrices by the use of Eq. ( l62l) . The results are depicted in Fig. [71 which 
shows the profile of perturbed azimuthal velocity in the rear and front of the planet, as well 
as the prediction of Eq. ( l62l) . We note that the jumps are not symmetric, although there is 
a good agreement between the theoretical prediction of the velocity jump on the one hand, 
and the average value of the rear and front velocity jumps on the other hand. One reason 
for that is that the perturbed pressure at the separatrices is not negligible compared to the 



0.40 
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Fig. 6. — Vorticity field in the vicinity of an Earth-mass planet at i = 64 orbits. The 
characteristics of the run are described in the text. The vorticity sheets are clearly apparent 
on the scparatriccs. Also, we note that the sign of the vorticity does not reverse at the 
stagnation point, but somewhere further on the downstream rear separatrix. The reason for 
this is that the perturbation of pressure is larger over the portion of the separatrix involved 
than at the stagnation point, as it crosses the wake. We also note that the vorticity is not 
singular on the upstream separatrices, which results from the fact that G is continuous across 
them. 
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Fig. 7. — Perturbed azimuthal velocity at t = 64 orbits, as a function of radius, at = 1 rad 
(dashed line) and at = — 1 rad (solid line). The vertical bar at the left shows the expected 
magnitude of the velocity jump, obtained from a measure of the torque excess using Eq. fl59|) . 
and using the result in Eq. (I62|) . 
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perturbed pressure at the stagnation point. Another one is the asymmetry of the horseshoe 
region: both separatrices do not have same value of 1(^2+ + Q_)/2 — Qp\. 

We note that for the run presented here, for which iS > 0, we have a negative torque 
excess, hence a negative pertur bation of surface density i n front of the planet, and a 



positive one behind the planet (IBaruteau fc Masset 



20081 ). The fact that we have a 



positive perturbation of azimuthal velocity in front of the planet (dashed curve of Fig. [7]) 
is compatible with an expansion of the material and consequently a decrease of surface 
density, because this perturbed velocity acts along with the Keplerian shear. Reciprocally, 
we have also a positive perturbation of velocity behind the planet, therefore opposed to the 
Keplerian shear, and compatible with a compression of material, hence an increase of the 
surface density. 

Finally, we mention that the vortensity is ill-defined at the separatrices. Whereas the 
vorticity is definable, the integral of its singular part yielding the azimuthal velocity jump, 
the vortensity would appear, on the separatrix, as the ratio of a singularity (that of the 
vorticity) by a non-continuous function (owing to the contact discontinuity, the surface 
density undergoes a jump at the separatrix), which is mathematically undefined. On the 
contrary, in the locally isothermal case (see paper I), the surface density field is continuous 
and the vortensity is defined everywhere. 



5.2. Production of vortensity 



On the contrary to a barotr opic situation, for w hich the vortensity is necessarily 



conserved along the streamlines ( iLovelace et al 



19991 ). some thermodynamical driving of 



the vortensity may occur during the horseshoe U-turns. We also note that far from these 
U-turns, the vortensity is conserved again, since all quantities depend only on r, hence the 
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pressure and density gradients are both radial, hence ahgned. 
Using Eq. fl2^ and (1271) . we are left with: 

drSo 



E 



UJ 



1-6T- 



(63) 



We comment that we have used the relationship dsG = —drSo/{4:ABx) to obtain Eq. (l63!l . 
The bracket of this equation corresponds indeed to material downstream of a horseshoe 
U-turn. During a U-turn, dcS is conserved, while x, to lowest order, reverses its sign, and 
Eq. fl23l) has been written for the unperturbed disk, that is to say upstream of horseshoe 
U-turns. The perturbation of temperature, within the horseshoe region, can be easily 
worked out assuming the conservation of entropy and an unperturbed pressure. This yields 



(iBaruteau fc Masset 



20081 ): 



6T 



-2T,S- 



(64) 



Using Eqs. (l63i) and (IMl) we get the production of vortensity during a U-turn: 









s 




UJ 


d 


UJ 



(65) 



2(7-l)A5r2 2B' 

We illustrate this in Fig. [HI which shows that there is a satisfactory agreement between 
the vortensity on the downstream side (here, outside, at the rear of the planet) and the 
vortensity predicted by Eq. fl65l) . We also note in this figure how the vortensity near 
corotation differs from the expected value, both on the inside and on the outside. On the 
outside, this is expected because the material very close to corotation did not have sufficient 
time to perform a horseshoe U-turn. On the inside, for 0.996 < r < 0.999, libration on 
tadpole streamlines brings vortensity from the outside, thereby yielding a perturbation of 
vortensity in the upstream region. More generally, the fact that streamlines are not exactly 
circular alters the initial profile of vortensity in this region. We also note the spike at the 



outer separatrix, resulting from the vorticity sheet presented in section 15. ![ where we also 
pointed out that, formally, the vortensity is ill-defined on the downstream separatrices. 
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Fig. 8. — Radial profile of S/u; at t = 48 orbits and = — 1 rad (solid line). The tilted 
dotted line shows the initial, unperturbed profile, while the tilted dashed line shows this 
profile symmetrized with respect to the corotation radius. It is therefore this profile that 
one would except if vortensity were conserved. If one takes into account the production of 
vortensity given by Eq. fl65|l . one obtains the dot-dashed line, which is in correct agreement 
with the observed profile. 
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We can make the following observations regarding Eq. fl65l) : 

• There is always a positive production of vortensity during the U-turns (hence a 
decrease of S/u;), regardless of the sign of S, and regardless of the location (at the 
rear or in front of the planet). 

• The production of vortensity is independent of the planet mass, and it occurs 
uniformly all over the downstream sides of the horseshoe region. While wider 
horseshoe U-turns bring the material closer to the planet, where the driving of 
vortensity is more efficient because the gradients of pressure and density are larger, 
they are performed faster than the more narrow U-turns, which occur further away 
from the planet. Eq. fl65|) shows that both effects cancel out, resulting in a fiat profile 
of vortensity perturbation. 

• If we translate the production of vortensity into a perturbation of density, assuming 
an unperturbed velocity field, then this perturbation has same sign on both sides 
of the planet. This is compatible with the fact that the bulk term of the horseshoe 
drag, i.e. the first term of Eq. (l37j) . is the same as in a barotropic situation. The 
weak thermodynamical driving of vortensity observed within the horseshoe region has 
therefore no impact on the torque. 

• The increase of vortensity after the first horseshoe U-turn is followed by a decrease of 
same magnitude at the following U-turn (in the idealized, inviscid situation considered 
here). This can be seen from Eq. fl27|) . in which we revert the d and u indexes. The 
vortensity of a given fiuid element therefore oscillates between its initial upstream 
value and its subsequent downstream value, and it does not undergo any drift on the 
long term. 
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5.3. Pressure disturbances 

In this section we wisli to determine the shape and amphtude of the pressure 
disturbances which arise in the coorbital region as the result of the horseshoe dynamics. 
In paper I, we have seen that the perturbation of vortensity was entirely accounted for 
by evanescent waves, so that the density or pressure response in the coorbital region was 
smoothed over a length-scale H . Here, we expect that the perturbation will be split into 
a localized part, which ends at the contact discontinuity at the separatrix, and a smooth 
part, corresponding to the launch of evanescent waves. 

Far from the planet, we assume purely circular motion, hence we can describe the 
disk's perturbed state with the perturbed azimuthal velocity 5v^{x), the perturbed surface 
density 5E(x), and the perturbed pressure 5P{x). We search the relationship between this 
set of perturbed functions, and the perturbations of vortensity 5w{x) and entropy 5s{x). 
We note that x has an arbitrary sign for the rest of this section. If we consider a radial 
profile in front of the planet, where the downstream separatrix lies inside of corotation, 
then Xs is meant to be a negative quantity. 

Using Eq. f[T21) . the disk's rotational equilibrium yields: 

-2Vl5v^ + + = Q (66) 

^0 ^0 ^0 

Anticipating on what follows, we note that while the radial derivative of an unperturbed 
quantity is of order ^o/'"; the radial derivative of a perturbed quantity 5^ is of order 

The second term of Eq. fl66|) is therefore r/H times larger than the third one. This 
equation can therefore be simplified as: 

d SP 

- 2no6v^ + ^ = (67) 
Denoting I = T^/u the inverse of the vortensity, we can write: 

61 = — - ^djv^, (68) 
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while the perturbation of entropy reads: 



5P 7 So 



5s =-^- -^5S. (69) 

^0 ^0 

Deriving Eq. fl67|) with respect to x, keeping only the derivatives of rapidly varying terms, 
i.e. the perturbations, and using Eqs. fIBSl) and (1^ so as to keep only the variable 5P, we 
obtain the following differential equation: 

5P_%.5P = c^So-, (70) 

where k, = (2f2oCc'o)^/^ is the epicyclic frequency, u = s^^'^l, and where we recall that 
is the adiabatic speed of sound (see section 12.11) . The general solution of Eq. (1701) is the 
convolution product of its right hand side by the Green's kernel K{x), which is the solution 
of: 

5P-%.5P = 5(x), (71) 

and whose expression is: 

K{x) = (27i/2/7)-ie-l^l/(^'''^), (72) 

where we have specialized to the Keplerian case, for which H = c^/n = c,/(7^/2^). This 
kernel represents the pressure response to a singular perturbation of w at a; = 0, of weight 
/ u{x)dx = Mo/(c^So)- We note that K{x) has a unitary weight: 

K{x)dx = 1. (73) 

Once the pressure response is known, one can infer the density response from Eq. fl69l) : 

. (If _ l£) . (r4) 

7 V^o So J 
Denoting 6T,d = SP/c^ and 6T,i = —T,o6s/{jso), we have: 



(5S = 6^d + S^i. (75) 
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In Eq. fl75l) . the subscript d stands for diffuse, since this component of the surface density 
scales with the perturbed pressure, which is smoothed over a length-scale H by the action 
of evanescent waves, as it involves the convolution of an arbitrary function by the kernel K 
defined at Eq. (1721) . Similarly, the subscript / stands for localized, since this component of 
the surface density scales with the perturbation of entropy, which is advected by the flow, 
and remains localized within the horseshoe region. 

Since 5P is the convolution product of K by c^So5u/u, we have: 

5Srf = ir*So— , (76) 

where * denotes the convolution product. Similarly, we have: 

ffi, = (77) 

The expression of the entropy perturbation is straightforward, as the latter is simply 
advected by the horseshoe dynamics. Since s oc r'^^ , we have 5s /sq = —2jSx/rp, hence 



(IBaruteau &: Masset 



20081): 

5Si = 2^oSx/rp. (78) 

The expression of the perturbation of u involves three components: the singular contribution 
of the vorticity sheet at the separatrix, the advection of the vortensity and entropy by the 
horseshoe dynamics, and the production of vortensity given by Eq. ( |65l) . We have mentioned 
in section 15.11 that the vortensity is ill-defined on the downstream separatrices, since the 
surface density is discontinuous there. Nevertheless, the jumps of surface density are of 
order Xg/a, hence to lowest order the vortensity can be estimated to be Av6{x — Xs)/So at 
the separatrices. We can therefore estimate the perturbation in u: 

(5 + V) + %^, (79) 



Su Av6{x — Xs) 2x^^ ^ ^ f^^prod 



Uq LOo Tp Iq 

where we have used the fact that u oc r"^"*"^, and where, using Eq. ( l65l) : 



lo 2{^-l)ABrf 
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For the sake of brevity, in what follows, we will keep the notation 5/prod, remembering that 
it corresponds to a uniform production of vortensity that scales with S"^, and which has 
same sign on both sides of the planet. The perturbation of density finally reads: 



5j: = k* 



-loAv6{x - Xs) So(5 + V) + 



Eq. (IHTj) features two kinds of terms: 



-So5. (81) 



• bulk terms (second to fourth term), which are defined for < x < Xg, prior to a 
possible convolution by K; 

• an edge term (the first one), which is defined at x = Xg, prior to the convolution by K. 

We better understand why the bulk of the horseshoe drag has the same expression as in 
the barotropic case. The convolution by K, which describes the spread of disturbances by 
evanescent waves, does not change the linear mass (J^ 6Ti{x)dx) of the perturbation, since K 
has a unit weight. Therefore, the torque exerted on the planet by the diffuse and localized 
stripes of perturbed density is the same as if all the disturbances remained localized, i.e. 
as if we omitted the convolution product in Eq. ( IHTl) . If we discard the singularity, which 
corresponds to the edge term (and, as we shall see, to the adiabatic torque excess), and 
the production of vortensity, which has same sign on both sides of the planet and does 
not affect the torque, we are left with the middle term of the bracket of Eq. (ISTIl . and the 
trailing, localized term, which simplify as: 

2r 

(5S = - — SoV, (82) 

exactly as in the barotropic case, prior to the convolution by the Green's kernel (see 
paper I). These properties are schematically depicted in Fig. O 

This first term of Eq. flHT]) corresponds to a single evanescent wave excited at the 
separatrix. It is the perturbed density associated with this wave that exerts the adiabatic 
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torque excess, since it is the only component that scales with the entropy gradient, owing 
to the simplification mentioned above. This is in agreement with Eq. fl57|) . which shows 
that the torque excess comes from an edge effect rather than a bulk effect, since it relies on 
a minute difference between the limit of the integration domain between the rear and front 
sides. We can check from Eq. ( l62ll that this edge term has correct sign. Since the perturbed 
density associated to this wave has a sign opposite of At>, the perturbed density has same 
sign as Fi in front of the planet, and an opposite sign behind the planet, as expected. 

A consequence of these properties is that, contrary to early expectations, it is 
not the localized component of the surface density that exerts the torque excess. The 
latter does exert a torque on the pl anet, which can be estima ted either by means of a 



partial horseshoe drag calcu^ 



(jPaardekooper fc Papaloizou 



at ions ( jBaruteau fc Masset 



20081 ). or by direct summation 



20081 ). It is found to have same sign as the excess, and same 
order of magnitude in usual setups. Nevertheless, this partial torque has to be added to the 
torque exerted by the diffuse component which scales as 5 + V, and the result is a corotation 
torque that scales only with the vortensity gradient. This can be expected on general 
grounds. The splitting of the perturbation into a localized component (an entropy wave) 
and pressure supported waves is very similar to the procedure used by numericists who use 
characteristic tracing to predict Riemann states, except that here the pressure supported 
waves are evanescent instead of propagative. The entropy perturbation regulates how much 
of the perturbation goes into the entropy wave, but it is the vortensity perturbation that 
dictates the total amount of mass perturbation, and it is therefore logical that the bulk of 
the corotation torque scales with the vortensity gradient. 

Owing to the simplicity of the functions involved in the convolution product, it is 
possible to write explicitly the expression of the pressure perturbation. This is done 
in appendix |Xl We show in Fig. [10] that the resulting expression, given by Eq. flASp . 
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Fig. 9. — Schematic representation of the different terms of Eq. (IHTI) . in their order of ap- 
pearance in the right hand side. The vertical arrow on the left plot represents the singularity 
at the separatrix. The first three terms represent the pressure supported density response, 
which is obtained by a convolution of the response depicted by the evanescent wave kernel 
of Eq. ( 1721) . The last term remains confined to the horseshoe region and is associated to the 
entropy wave. The terms are here depicted for x > 0, i.e. behind the planet. In front of 
the planet, they would have an opposite sign, except the third term, corresponding to the 
production of vortensity, which would remain negative and has therefore no impact on the 
torque. 




Fig. 10. — Profile of pressure perturbation for the calculations EAq (a), EAg (6), EAig (c), 
EA24 (d), EA32 (e) and EA40 (/). The solid line shows the profile observed in numerical 
simulations, evaluated using Eq. ( 183|) . and the dashed line shows the result of Eq. ( lASp . 
The dotted line shows the partial result of Eq. (lASp . in which we omit the first term, which 
corresponds to the isolated evanescent wave at the separatrix corresponding to the torque 
excess. The results of the plot d have been multiplied by —3, in order to make it more legible. 
This case corresponds to an approximate cancellation of the pressure term stemming from 
the S + V term and of the single evanescent wave excited at the separatrix. The case of the 
plot e is also of interest: for this calculation, V + iS ~ 0, as can be seen from the dotted line, 
hence the pressure perturbation almost exclusively consists of the single evanescent wave 
excited at the separatrix, which is responsible for the torque excess. 



satisfactorily reproduces the perturbations of pressure observed in numerical simulations. 
Namely, we transform the profile of perturbed pressure from a numerical simulation as 
follows, in order to filter out the term arising from the production of vortensity, which is 
poorly reproduced by a rectangular function, as it has been shown in Fig. [H) 



We then evaluate the theoretical expression of the perturbed pressure given by Eq. (lASp . in 
which we use for Fi the estimate provided by Eq. fl59|) (corresponding to the thick dashed 
curve of Fig. [S]), and for Xg the arithmetic mean of Xp and x^, determined by a streamline 
analysis. Owing to the symmetrization performed in Eq. fl5^ . we discard the term in 5/pi.od 
in Eq. (1A5I) . The correct agreement between the observed and predicted pressure profile 
indirectly confirms that it is indeed the single evanescent wave excited at the separatrix 
that is responsible for the torque excess. We also recover a feature noticeable in Fig. O 
While the plots a and / of Fig. [10] correspond to the same absolute value of the entropy 
gradient (namely 8/7), the single evanescent wave has a much larger amplitude in the case 
iS > (i.e. for the plot /), almost of a factor of two, as can be seen from the distance 
between the dotted and dashed curves. This is compatible with the fact that, in Fig. [5l we 
found that the corrected torque excess, in absolute value, is almost a factor of two larger at 
S = 8/7 than at 5 = -8/7. 



We have seen in paper I how the excitation of evanescent pressure waves alters the 
width of the horseshoe region, in a different manner in front of the planet and behind the 
planet, so that it renders the horseshoe region asymmetric. This asymmetry, however, 
was found to have virtually no impact on the horseshoe drag, while the separatrices all 
shared the same value of the Bernoulli invariant. In the adiabatic case considered here. 




(83) 



6. Interpretation: an intrinsic asymmetry 
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the evanescent waves launched downstream of the U-turns (see section 15.31) also alter the 
horseshoe width and render the horseshoe region asymmetric, but there is another source 
of asymmetry that preexists the launch of evanescent waves. Let us assume for a moment 
that the disk is composed of non-interacting test particles, so that there are no evanescent 
waves excited by the horseshoe motion, but let us also assume that the separatrices of the 
horseshoe region are given by the G+ and G_ values, as in an adiabatic disk (the expansion 
of G to lowest order in x, given by Eq. f|T9|) . is independent of the disk's thermodynamics.) 
A different value for G_ and implies a different upstream half width of the horseshoe 
region. Such a case is depicted in Fig. [HI The test particle A, which lies outside of the 
separatrix, is circulating, and therefore mapped to A' by the flow. Similarly, the test 
particle B lies inside of the separatrix, and thus is librating. It is therefore mapped to B' 
by the horseshoe motion. There is therefore an over-density in the darker stripe in the 
bottom right quadrant of Fig. [TT|, as test particles of different origins merge into this stripe. 
Reciprocally, a void region appears at the downstream rear separatrix (the white stripe in 
the top left quadrant of Fig. [TTi) . This region can neither be reached by the test particles 
executing horseshoe U-turns, as these are too narrow, nor can it be reached by test particles 
of the outer disk. The net effect of the asymmetry is that two stripes of perturbed surface 
density iSg appear at the downstream separatrices, with a width 



5x., 



AG 




AG 


drG 




AABxs 



The linear mass of these stripes is therefore, in absolute value: 

AG 



A 



AABx. 



^|At;| = /o|At;|, (85) 



2B 

where we have used Eq. (1621) . This linear mass corresponds to the factor of the 5-function 
in the bracket of Eq. (ISTI) . The torque due to the stripes can be evaluated as follows. The 
symmetric, barotropic case would be recovered if one sent the excess of surface density of 
the overdense stripe to the empty stripe, spending angular momentum for this purpose at 
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the rate: 

iFil = AaBxs ■ \2Axs^o6xs\, (86) 

where the factor in absolute value represents the mass flow rate in the stripes. Using 
Eqs. f l8^ and fl86|) . one recovers the expression of the torque excess, given by the last term 
of the right hand side of Eq. fl37j) . In the case considered in Fig. [11], which has a wider 
horseshoe region at the front of the planet, one would have to give angular momentum 
to the material in excess in the front stripe, in order to recover a symmetric situation. 
Differently said, this implies that in the adiabatic situation the disk's material receives less 
angular momentum than in the barotropic case, hence that there is a negative torque excess 
on the disk, and a positive torque excess on the planet, by virtue of the action and reaction 
law. This can also be deduced simply from the distribution of perturbed density, since there 
is a positive perturbed density in front of the planet and a negative one behind the planet. 

In a disk with pressure, narrow stripes of surface density are spread radially by the 
excitation of evanescent waves, which are not easily detectable (see section [573]) . They leave 
however an imprint on the flow in the form of vorticity sheets. These are the most tangible 
perturbations of the flow associated to the torque excess. 

One can also understand that the torque excess due to these stripes can in principle 
be extremely high. While the surface density perturbation in the barotropic case typically 
amounts to 0(T,QXs/a) (prior to the convolution by the evanescent wave kernel), the stripes 
contemplated here correspond to a perturbation of surface density that is iSg. If the 
asymmetry becomes large enough to represent a sizable fraction of Xg, the linear mass of the 
perturbation associated to the excess supersedes the linear mass of the vortensity related 
perturbation by a factor 0{a/xs), that is to say typically by two orders of magnitude, for 
the situations considered here. 

Further insight into the origin of the asymmetry can be gained by considering, for the 
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> 

A zimuth 

Fig. 11. — Schematic representation of an asymmetric horseshoe region. The dotted lines 
represent the symmetric separatrices of a barotropic situation. With respect to this case, 
the rear separatrix is shifted towards corotation, while the front separatrix is shifted away 
from corotation. The width Sxg of the stripes has been exaggerated to improve legibihty. 
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sake of simplicity, an idealized situation in which the stagnation point lies at = 0, i.e. at 
the maximum of the perturbed pressure. We can then evaluate, in order of magnitude, how 
the torque acting on a fluid element moving along an upstream separatrix differs from the 
barotropic case, and infer from this how the separatrix should be shifted. This situation is 
depicted in Fig. [12l For the sake of definiteness we assume that there is a negative radial 
entropy gradient. 



see 



The perturbation o: 



Baruteau fc Masset 



density associated to the advection of entropy in the adiabatic case 



20081 ) enhances or diminishes the specific pressure torque — VP/S. 
As can be seen in Fig. [T21 a consequence of this change is to render the horseshoe region 
asymmetric: the rear part shrinks, whereas the front part expands. A coarse estimate of the 
variation of the width of the horseshoe region can be given as follows. A fluid element on 
a separatrix, when it reaches the azimuth of the stagnation point, has received exactly the 
amount J of angular momentum necessary to bring it to the stagnation point (by definition 
of the separatrix): 

J= (-d^<l>-^]dt = 2BrpXs, (87) 



where ts is the date at which the fluid element reaches the stagnation point. Assuming 
the fluid element to have at any instant in time the azimuthal velocity of the unperturbed 
Keplerian shear, we can write: 

A rough estimate of how Xs varies when the pressure term is modified can be obtained 
by taking x out the integral in Eq. (188|) . assuming that it is everywhere equal to Xs- 
At this level of approximation, the perturbation of the pressure term can be written as 
~ (90P/So)(25x,/a). We are then left with: 



i^f' .(25a;,/a) ~ 2ABxs5xs, (89) 
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Fig. 12. — Schematic view of the action of the pressure gradient in an ideahzed case. The 
grey contours represent the perturbed pressure. The sohd hues depict the separatrices of a 
barotropic case. The grey thin arrows represent — VP/E in the barotropic case, and the black 
thin arrows the same quantity in the adiabatic case, the difference arising essentiaUy from 
the variations of E between both cases. The thick arrows represent the excess of — VP/S, 
which is negative on both sides of the planet. The fluid particle initially located on the rear 
upstream separatrix of the barotropic case therefore becomes circulating, as depicted by 
the dashed streamline. Similarly, the fluid parcel of the front separatrix loses more angular 
momentum than it would lose on the barotropic separatrix, hence the streamline associated 
to it is librating, as indicated by the dotted streamline. 
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which can be recast as: 



Sx, ^ . (90) 



We note that we obtain the same expression for the variation of the horseshoe width using 
Eqs. fHTl) . fl50l) . fl54l) and fl8^ . within a factor of order unity (the value 5xs obtained above 
just refers to the variation of the width of one side of the horseshoe region, whereas the 
expression worked out earher in this section and depicted by the thin double arrows in 
Fig. [11] corresponds to adding the width variation of the front and rear sides.) 

This order of magnitude estimate allows to identify the origin of the horseshoe 
asymmetry, and allows to understand why the adiabatic torque excess scales with the 
perturbation of pressure at the stagnation point (p^), at least in the regime of small 
planetary masses (for which ps <^ Pq.) We finally note that, while this order of magnitude 
was obtained in the simplified case of a stagnation point located at the maximum of 
perturbed pressure, it still holds when the stagnation point is azimuthally shifted, as in 
real situations. A more refined treatment of the dynamics of a fiuid element moving along 
the separatrices would have exhibited its Bernoulli invariant, and it would have led to an 
exact relationship such as Eq. fj44|) . which is equivalent to Eq. flHHl) . This relationship shows 
that the width of the horseshoe region depends exclusively of the fiow properties of the 
stagnation point, independently of the location of the latter and of the path followed by the 
fiuid elements. 



7. A suitable torque expression 



The torque expression of Eq. (jSSD involves the planetary potential at the stagnation 
point, as well as the distance of the upstream separatrices to corotation. As such, it is not 
well suited for torque estimates, and needs further transformation. 
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We have seen in section 14.21 that the half width of the horseshoe region has a complex 
dependence on the disk's and planet's parameters in the adiabatic case. In particular, at 
large values of \S\, the horseshoe region is wider than in the barotropic case, which boosts 
the vortensity related corotation torque. The global torque excess is therefore the sum of 
two non-trivial individual excesses: the one related to the gradient of entropy, which has 
been studied in depth in the preceding sections, and the one related to the boost of the 
vortensity related corotation torque. We note however in Fig. [5] that the global excess, 
which corresponds to the dotted curve, has an almost linear dependence on the entropy 
gradient. Our aim is therefore to exhibit a separable form of the torque excess that involves 
the product of iS by a function of e. This expression aims at describing the global torque 
excess, therefore accounting for the whole difference between the adiabatic and barotropic 
case. For this purpose we can simply use the data of the series SAj, which has a null 
vortensity gradient and for which the adiabatic torque excess accounts for the whole excess. 
We then rescale the function obtained to get the excess at any entropy gradient. 



7.1. Scaling of the torque excess 

We first establish the scaling of the torque excess as a function of the disk's and planet's 
parameters. Eq. fl55l) shows that the torque excess involves the half width of the horseshoe 
region, and the potential at the stagnation point. In order to remain in the framework of a 
separable expression, we neglect the dependence of Xg on the entropy gradient, mentioned 
in section 14.21 and we use the barotropic approximation (i.e. for 5 = 0) to infer the half 
width. Using Eq. (157|) . we obtain: 

- tI"/!- pi) 

where C(i) is a dimensionless coii stant that depends on the softening length 



( jPaardekooper fc Papaloizou 



2009bl ). Examination of Fig. [3] shows that this approx- 



imation leads to errors of at most ~ 30 %. We check in what follows that the potential 
at the stagnation point can also be considered as essentially independent of the entropy 
gradient. 



7.1.1. Potential at the stagnation point 

Evaluating the planetary potential at the stagnation point requires to know the location 
of this point. Writing the Euler equations in steady state, in a system of units in which the 
unit of length is H = Cg/Q, one can realize that, as long as the planet is deeply embedded 
(i.e. aq^^^ <^ H), the azimuth of the stagnation point (which lies almost on corotation) has 
to scale as hf{e), where / is a function that we shall determine. We have therefore: 

= (92) 

Fig. [13] shows the location of the stagnation point as a function of e in an adiabatic flow 
with 7 = 1.4, obtained from the series of numerical simulations EAo<j<2o (see section UT]), 
and a similar series with a different entropy gradient. Under e ~ 1.0, the distance of the 
stagnation point to the planet is larger than the softening length, and the denominator 
of Eq. (!92|) is dominated by the function /(e), which is found to hardly depend on S, as 
both series essentially coincide for e < 1.0. Reciprocally, for e > 1.0, the stagnation point 
lies within a softening length from the planet, and the potential at the stagnation point 
approximately scales as e~^, and is therefore also independent of the entropy gradient. 

7.1.2. General expression of the torque excess 
Using Eqs. , and , we obtain: 

Ti = ~A{e)STo, (93) 
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Fig. 13. — Location of the stagnation point as a function of the softening length, measured 
at t = 40 orbits. Diamonds show values obtained in the EAo<j<2o series. The two squares 
show values obtained for the same parameters as in the EA series, except that the resolution 
has been increased to 2000 x 2000. The stars show values obtained from a series similar to 
EA, except that it has a different entropy gradient: S — +12/7. The dotted curve, which 
shows [/(e)^ + e^]^/^, has been obtained using the data provided by the EA series. 
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where Tq = TioVtiq^ / h? , and where: 



A{e) 



ACHt) 



1.4 



(94) 



where we have speciahzed to the Keplerian case {A = —3Q/A, B = Q/4). Note that we 
exphcitly restrict ourselves to the case 7 = 1.4. Ahhough Eq. (19T|) apphes for arbitrary 
values of 7, the azimuth of the stagnation point, considered in section [7.1.11 may depend 
on the adiabatic index. Keeping an explicit dependence in 7 in Eq. (l94l) could therefore 
be misleading. The adiabatic torque excess, given by Eq. (1^ . h as a scaling similar to the 
barotropic corotation torque or the differential Lindblad torque ( Tanaka et al.|[2002l ). The 
only differences with the barotropic corotation torque is that the excess scales with S rather 
than V, and that it involves a different dimensionless function of e. 

In order to determine this function, we use the series SAj, as explained above. We 
restrict ourselves to the range 0.1 < e < 1. Values of e larger than one are not relevant, 
whereas values smaller than 0.1 are difficult to study owing to the very high resolution 
required. As mentioned in section I^^Sj the torque excess exhibits a dependence in with a 
good approximation over the range of e considered (a dependence partly due to the increase 
of Xs at small softening length, and partly due to the increase of the planetary potential at 
the stagnation point at small softening length). Fitting the torque excess over the range of 
e considered, we obtain: 

A{i) ^ ^. (95) 

The expression of the torque excess is therefore: 

1.35, 



(96) 



One can see in Figs. H] and O that this expression (which is represented by the grey dashed 
line) is in good ag reement with the global ad iabatic excess. It is also in good agreement 



with the results of 



Baruteau fc Masset 



fl2008f ) 
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7.2. Considerations on the three-dimensional case 



If one considers that a two-dimensional case with softe ning length e represents the 



Masset 



2OO2I ). one can perform 



layer at altitude z = ±e of a three dimensional case (see e.g. 
an integration of the torque over e in order to get an approximate value of the three 
dimensional torque. A layer comprised between altitude z and z \ dz has a surface density 
(So/\/27rif) exp[— 2;^/(2if^)](i2;, hence the three dimensional adiabatic torque excess should 
be approximately given by: 

-00 V' 



3D 



exp[-z7(2i7^)]A . ^ 



(I) rldz 



-\ -SY. 



A{t)de 



(97) 



While this kind of calculation can be performed relatively easily for the vortensity 
related torque (for which the dimensionless function of e converges when e ^ 0) or for the 
differential Lindblad torque, it is a risky exercise in the case of the adiabatic torque excess, 
since the approximate expression of Eq. ( l95l) diverges for e — 0. For the setup that we 
consider, we have tried to perform calculations at higher resolution and smaller smoothing 
lengths. The smallest smoothing length for which we have a steady situation is e = 0.07i/ 
(corresponding to the left square in Fig. [T^ . Under this value, we end up with very messy, 
unsteady situations, in which the flow is strongly affected by the presence of numerous 
vortices which drift along the separatrices. We can therefore only provide a conservative 
estimate of Eq. (197|) . by truncating the integral at e = 0.07. We obtain: 

Vf^ ^ -2.85ro. (98) 

In order to estimate the total torque in the three-dimensional case, we need to know the 
differential Lindbla d torque, and the vor tensity related part of the horseshoe drag. The first 



of those is given by 



Tanaka et al 



(I2OO2I ). and is, with our notation: 
ri^ = -(2.34-0.1a)ro. 



(99) 
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for a disk with out a temperature K radient. We note that we do not use the total torque 



value given by 



Tanaka et al. 



(I2OO2I ). but only the differential Lindblad torque. The total 



torque value contains indeed t he linear corotation torque, which i s not the corotation 



torque exerted in steady state (jPaardekooper fc Papaloizou 



2009al ). Instead, we consider 



separately the vortensity related horseshoe drag, and apply to it a treatment similar to 
that of Eq. 0971) . which is straightforward in this case, because of the absence of divergence 



when e ^ 0. The two-dir aensional component of t 



vortensity gradient reads (IWard 



1991 



Masset 



le horseshoe drag that scales with the 



2OO1I): 



= -VC\e)T,, 



hence its three-dimensional estimate is: 

3 



p3D 
^ V 



57r 



Jo 



(100) 



(101) 



where we have used the series Slj to tabulate C(e). The total torque acting on a low-mass 
planet in a three-dimensional disk reads therefore: 



r = [-(2.34 - 0.1a) + 0.93V - 2.85]ro, 



(102) 



where we recall that the last term is a conservative estimate, and is likely to actually have a 
numerical coefficient larger than 2.8. We see that for |iS| ~ |V| ~ 1, the term associated to 
the entropy gradient dominates the horseshoe drag, hence the adiabatic torque excess should 
efficiently halt or reverse mi gration in disks with mode r ate, ri egative entropy gradients, in 



agreement with the result of 



Paardekooper fc Mellema] (120061 ) . who found that planetary 



migration can indeed be reversed in a three-dimensional disk. 



8. Discussion 

In this section we describe some side results, and we draw a list of points that deserve 
further investigation. 



8.1. Topology of the flow and location of the stagnation point 



As stated in section \2A\ we have assumed in this work that there is only one X-point 
in the vicinity of the planet. We find that this is indeed the case in most situations. We can 
nevertheless get two X-points in the vicinity of the planet, if one of the following conditions 
is fulfilled: 

• the entropy gradient vanishes, or is very small, 

• the softening length of the potential is very small, 

• the fiow has not reached a steady state and is observed before a U-turn time scale. 

While the first and last conditions are not relevant to the present work, which assumes a 
steady state and a sizable entropy gradient, the second case has some importance, as it 
gives indications of what can happen near the equatorial plane in a three-dimensional case. 
In the cases in which we have two X-points that subsist during the whole simulation, we 
get very time-dependent torque values. This can be interpreted as entropy trapped in the 
small libration island defined by the two X-points. The perturbation of density associated 
to the entropy distribution, that librates in this island, gives rise to large torque variation, 
as it corresponds to material very close to the planet. We note that this is quite in contrast 
with what we have claimed in section 15. 3[ where we have argued that it is the vortensity 
perturbation that yields the bulk torque. Whereas this is true for the downstream horseshoe 
stripes which lie at some distance of the planet, here the planet lies within the libration 
island, and the fact that the density response is spread radially or not definitely matters. 
We also note that the present discussion offers a number of similarities with the discussion 
in paper I about the fiow topology as a function of the temperature gradient and of the 
softening length, except that here it is the entropy gradient that plays a role, rather than 
the temperature gradient. 
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When the flow is observed at an early stage, it may display two X-stagnation points, 
which are generally not on the same separatrix. We observe that once the flow has settled 
in a steady state, it is the outermost stagnation point (the one that is on the separatrix 
that lies further from corotation) that subsists. It turns out that this point lies behind the 
planet for S < 0, and is in front of the planet for S > 0. While this fact is unimportant 
for a planet in flxed circular orbit, where only the distance of the stagnation point to the 
planet matters to evaluate the torque excess, it is crucial when the planet migrates, as 
the stagnation point is shifted from the position that it has on a flxed circular orbit. The 
consequences of this are discussed in the next section. 



8.2. Feedback on migration 

Since the adiabatic torque excess depends on the flow properties at the stagnation 
point of the horseshoe region, and since, as we shall see below, the location of the latter 
depends on the drift rate, we expect that a feed back loop can be e stablished between the 



total torque and the migration rate, much as in type III migration (IMasset Sz Papaloizou 



20031 ) . Contrary to type III migration, however, the feed back should generally be negative, 
therefore acting to decrease the absolute value of the drift rate. It should also be of 
relatively minor importance, except in very massive disks, with Toomre's Q parameters 
close to unity. Consider the following numerical experiment: we repeat the runs SAio and 
SIio, except that we impose an inwards disk drift = — 1.9 ■ 10~^. We note that, if the 
planet was allowed to freely migrate in run SAio, it would do so at the rate d = 1.9 • 10"^. 
Rather than releasing the planet, we impose an inwards drift of t he disk, of same magnitude . 
This technique has already been used for type III migration (see lMasset fc Papaloizoull2003l . 
section 5.6). We see in Fig. [H] that the torque in the isothermal case does not depend on the 
drift rate, while there is a noticeable difference in the adiabatic case, and that the torque 
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tends towards a somehow smaller value. The reason for this can be understood from the 
profile of radial velocity at corotation, depicted in Fig. [T51 If the planet migrates outwards, 
the stagnation point (in a frame that moves radially with the planet) is located where the 
horseshoe dynamics endows fluid elements with a radial velocity equal to that of the planet 
(similarly, if the planet is kept fixed, the new stagnation point is expected to lie where the 
radial velocity in the initial case is — r^.) Since fluid elements move outwards behind the 
stagnation point, the azimuth of the stagnation point decreases, hence the latter recedes 
from the planet, since it is located behind the planet in the case considered. A consequence 
of this recession is a lower perturbed pressure at the new position of the stagnation point, 
hence a lower value of the adiabatic torque excess. A similar decrease (in absolute value) 
is also expected in the case S > 0. The stagnation point is then in the front of the planet 
(see section I5TT]) . and the total torque acting on the latter is negative. The stagnation point 
in the migrating case is therefore expected where the horseshoe U-turns are performed 
inwards, that is to say at the front of the stagnation point of the fixed case, and therefore 
further again from the planet. We have checked this with additional simulations, not 
reproduced here. 

A number of comments are in order: 

• Fig. [15] shows a poor agreement between the radial velocity for the disk's drift case 
and the original radial velocity offset by — r^, in the region of the stagnation point 
(—0.07 < (p < —0.03). This is presumably due, at least in part, to the resolution of 
the grid. The direction of the shift of the stagnation point, and its order of magnitude, 
are nevertheless compatible with the new stagnation point being determined by the 
location where, in the initial flow, Vr = a. 

• Since both Vr and a scale with the planet mass (for the range of planetary masses 
considered here), the shift of the stagnation point should not depend on the planetary 
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mass. 

• Since a scales with the disk mass, the shift of the stagnation point of a freely migrating 
planet should increase with the disk mass. The magnitude of this effect, for the disk 
considered in the numerical experiment described above, is at about 1/20*^ of the 
total torque. The disk has a Toomre parameter Q ~ 8. This shows that the effect we 
describe is of marginal importance, except in the most massive disks, for which Q ^ 1- 

• The drift rate in steady state is given by: 

2Sad = 7+^d, (103) 
oa 

where we assume that the torque excess, to lowest order, has an affine dependence on 
d, and where 7 is the specific torque acting on the planet. We therefore obtain: 

In the particular case described above, we have d^/dd ^ —0.05 • 2Ba, hence the 
steady state drift is almost equal to that dictated by the torque value measured in the 
disk drift case (i.e. further iterations with a drift rate of the disk given by the new 
torque value would hardly change the result). 



8.3. Extension to arbitrary temperature profiles 

The analysis presented in this work suffers from a restriction very similar to the 
restriction of the isothermal case, which is that the temperature profile is assumed to be 
fiat. We have nevertheless performed many additional calculations in which we relaxed 
this constraint, and found that our main results essentially hold in disks with arbitrary 
temperature profiles: 
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Fig. 14. — Specific torque on planet in the runs with an inwards disk drift. The solid line 
shows the adiabatic case and the dotted line shows the isothermal case. Grey curves show 
the initial runs (without disk drift). The inset plot shows a close up on the curves (the value 
of the torque at t = 21 orbits in the case without disk drift has been subtracted to improve 
legibility) . 



-65- 




-2.0x10'^ I ^ ^ ^ \ ^ ^ ^ L ^ ^ ^ ^ ^ ^ \ ^ ^ ^ 

-0. -0.08 -0.06 -0.04 -0.02 0.00 

A zimuth 



Fig. 15. — Radial velocity as a function of the azimuth, at r = Tc, in the static case (dashed 
line), and in the case with disk drift (solid line). The grey dashed curve shows the curve of 
the static case offset of fa- The rightmost vertical dotted line shows the initial location of 
the stagnation point, while the leftmost one shows the location of the stagnation point in 
the case with a disk drift, which is also the location expected for a planet moving outwards 
in the non drifting disk. 
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• The global torque excess, defined by Eq. fl571) . scales with the entropy gradient. 
The isothermal runs which we use to apply Eq. flFTI) in the general case have same 
profiles as their adiabatic counterparts, hence they are locally isothermal runs, which 
have singular evanescent waves excited at the separatrices, contrary to the globally 
isothermal runs (see paper I). These components, which scale with the temperature 
gradient, should therefore be also present in the general, adiabatic case. 

• The half width of the horseshoe region exhibits the same characteristic V-shaped 
dependence on S, as shown in Fig. [3l 

The adiabatic torque excess, as a function of the entropy gradient, displays little scatter, 
and is therefore essentially a one-to-one correspondence. The results presented here can 
therefore be used to predict the torque in a general situation. They indicate that, in the 
general case, the horseshoe drag consists of three terms: the classical vortensity related 
horseshoe drag, the adiabatic torque excess presented in this work, and an additional term 
which scales with the temperature gradient (see paper I). The latter corresponds to the 
creation of vortensity in the vicinity of the planet, near the stagnation point, where the 
fiow is slow and the driving of vortensity is efficient. We comment that the adiabatic torque 
excess can also be seen as a singular creation of vortensity at the downstream separatrices, 
with the word of caution that vortensity is mathematically ill-defined at these separatrices, 
because of the presence of a contact discontinuity (so that one has to consider the low order 
estimate cu/Sq, rather than cu/S). 

8.4. Adiabatic excess in turbulent disks 

The analysis presented above relies upon the existence of a steady state and on that 
of a unique stagnation point. It therefore does not apply to the case of a turbulent disk. 
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The adiabatic excess has been interpreted, however, as due to a constitutive asymmetry 
of the horsehoe region, acquired by the upstream flow under the action of the pressure 
gradient on overdense or underdense regions. To some extent, this effect should therefore 
persist in a turbulent disk, since we expect to find, in a turbulent disk, regions that are, in 
average, under- or over-dense with respect to the barotropic case, as a consequence of the 
advection of entropy. This issue deserves significant further work, since most of the regions 
of the disks where planetary migration is supposed to take place should be turbulent. In 
particular, it would be of interest to know whether the time-averaged torque excess is equal 
to the steady state estimate. 



8.5. Saturation issues 

In the absence of any process that can transfer angular momentum between the 
horseshoe region and the rest of the disk, the horseshoe drag is bound to saturate. A 
non-vanishing time averaged value of this torque would indeed imply a sustained transfer 
of angular momentum from the planet to the horseshoe region, where it would accumulate. 
Viscous torques, exerted at the separatrices between the ho rseshoe region and the outer or 



20011 1. In the case of 



inner disk, can ensure this transfer of angular momentum (iMassetl 
an isothermal disk, viscous diffusion tends to restore the large scale gradient of vortensity 
across the horseshoe region. If it can do so in less than a libration time, the vortensity 
upstream of the horseshoe U-turn is that of the unperturbed disk, and the horseshoe drag is 
permanently equal to its initial value. In the case of an adiabatic flow, the vortensity related 
part of the horseshoe drag follows the same picture, and a sufficiently large viscosity is 
required to avoid the saturation of this component of the drag. The adiabatic torque excess, 
however, results from the advection of entropy, and to the existence, at the stagnation point, 
of a discontinuity of the entropy brought from the inside and from the outside. Sustaining 
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the torque excess th erefore implies that the discontinui ty of entropy at the stagnation 



point is maintained. 



Paardekooper fc Papafoizou 



(I20081) consider therm al diffusion as the 



Kley fc Cridal (120081 ) consider realistic 



dissipative process that fulfills this function, while 
radiative effects. In both cases there can be a sustained, positive torque acting on the 
planet. We note however: 



• that in any case a finite amount of viscosity is required. Thermal diffusion or radiative 
effects, which do not feature in Euler's equation, cannot ensure the transfer of angular 
momentum out of the horsehoe region. 

• That there is no theoretical expression of the steady state horseshoe drag that takes 
into account the balance between saturation and the dissipative processes that prevent 
it. Such an expression is definitely required for planet population synthesis based on 
one-dimensional models of migration. 



We also comment that, while saturation can be interpreted as the phase mix- 
i ng of the contribution of diffe rent horseshoe streamlines in the barotropic case 



(IBalmforth fc KorycanskyI 120011 ). the case of the adiabatic torque excess is slightly different 
because there is essentially one streamline involved in this excess: the separatrix. It is 
therefore possible that the torque excess oscillates over a very long time, or indefinitely. 
Numer ical simulations show that the oscillations of the torque damp after a few libration 



times ( iBaruteau 



20081 ). It is not clear, however, whether this is a physical effect or the 
consequence of numerical diffusivity. This again appeals for the need of simulations at very 
high resolution. 

We finally mention the systematic appearance of a vortex at one of the downstream 
separatrices (the one that has a negative singularity of vorticity). It would be of interest to 
study the impact of this vortex on the saturation of the adiabatic excess. On the one hand. 
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it provides a mixing of the entropy in the vicinity of the separatrix, and on the other hand, 
it may assist the horseshoe region in exchanging angular momentum with the rest of the 
disk. 

9. Conclusions 

We have shown that the horseshoe drag exerted on a low-mass planet, embedded in a 
disk that behaves adiabatically on the time scale of a horseshoe U-turn, can be considered 
as the sum of two terms: 

• a bulk term, involving all the horseshoe streamlines, which necessarily scales with 
the vortensity gradient, exactly as in a barotropic disk. This is in agreement with 
the forecast made in the introduction: there is no room, for a generic horseshoe 
streamline, for an exchange of angular momentum with the planet different from that 
of a barotropic situation. 

• An edge term, which manifests itself as a torque excess with respect to the barotropic 
case. This term scales with the radial gradient of entropy, and it corresponds to a 
constitutive asymmetry of the horseshoe region. In a pressureless disk, this asymmetry 
would trigger the appearance of underdense or overdense thin stripes of gas at the 
downstream separatrices. In a gaseous disk, pressure forces spread radially these 
stripes by the excitation of evanescent waves, leaving imprints as vorticity sheets at 
the separatrices as the most tangible manifestation of this effect. 

We provide an expression for the torque excess, given by Eq. ( |96ll . which should be added 
to the differential Lindblad torque and to the vortensity related horseshoe drag whenever 
the flow behaves adiabatically on the timescale of the horseshoe U-turns. We also provide a 
tentative expression of the total torque in the three-dimensional case, given by Eq. fll02p . 
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We find that the origin of the adiabatic torque excess is not the conspicuous under- 
or over-dense regions that appear within the horseshoe region and which are bound by a 
contact discontinuity at the downstream separatrices, as originally thought. Instead it is 
due to rather discreet single evanescent waves launched at the downstream separatrices. 



This answers some issues arising from earli er works. In particu 



excess is expected in the limit of a cold disk (iBaruteau fc Masset 



ar, n o adiabatic torque 



20081 ). whereas the contact 



discontinuities still exist in this limit and delineate perturbations of finite mass. The torque 
excess that we find scales with the perturbed pressure at the horseshoe's stagnation point, 
hence it vanishes in the limit of a cold disk. 

We found a number of side results: 



• The horseshoe region in the adiabatic case has a more complex dependence on the 
disk's parameters than in the barotropic case. In particular, the horseshoe region is 
wider than expected from barotropic estimates, when there is a non-vanishing entropy 
gradient. A consequence of this effect is that the bulk horseshoe term, that scales with 
the vortensity gradient, is boosted with respect to its value in a barotropic disk. This 
effect is relatively weak, however, and most of the difference between the adiabatic 
and barotropic case comes from the edge term, in disks with realistic profiles. 

• The adiabatic torque excess slightly depends on the migration rate. There is therefore 
a weak feed back on migration. This feed back is found to be generally negative (i.e. 
it tends to lower the drift rate, either inward or outward), and is virtually negligible 
except in very massive disks. 



We find that even at relatively high resolution, numerical simulations may mistakingly 
locate the stagnation point, which has some impact on the torque excess value. We therefore 
stress the need for very high resolution calculations, for which nested mesh codes would 
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be a valuable tool. In a similar vein, the steep dependence of the excess on the potential's 
softening length suggests that the effect can be very strong in the three dimensional case. 
This issue requires significant further work. 

The numerical simulations performed in this work have been run on the 92 core cluster 
funded by the program Origine des Planetes et de la Vie of the French Institut National 
des Sciences de I'Univers. Partial support from the COAST project {Computational 
ASTrophysics) of the CEA is also acknowledged. The authors also wish to thank G. 
Koenigsberger for hospitahty at the Instituto de Ciencias Fisicas of UNAM, Mexico, and 
acknowledge partial support from CONACYT project number 24936. The authors are 
grateful to S. Fromang and C. Baruteau for a thorough reading of a first draft of this 
manuscript. 



A. Expression of the pressure perturbation 



The convolution product of the rectangular function 



R : X t-^ 1 ii < X < Xo 



(Al) 



otherwise. 



by the kernel K{x) — exp( 



a;|/A)/(2A) is: 




(A2) 




otherwise. 



Similarly, the convolution product of the triangular function 



T:a;i— > a;ifO<a;<a;o 



(A3) 



otherwise. 
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by the kernel K{x) = exp(— |a;|/A)/(2A) is: 

exp(x/A) 



T(x,xo, A) 



xq + A) exp ( — J-] + A 



if X < 0, 



(A4) 



exp (—a;/ A) 



[xo - A) exp ( y ) + A 



x + -exp 



x\ Xo + A 



exp 



if X > Xo 

X — Xo 



A 



if X G [0, Xo]. 



Using Eq. (IHTj) . we obtain: 



6P 



exp 



(8a|A|Sx2). (271/2/7) 

+ SoC^^^y^^(x, Xs, y/lH), 



X 



2Soc2 



{S + V)f{x,Xs,^H) (A5) 



where we have also made use of Eq. ( 1621) to transform the first term. 
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